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We report the dependence of Ruderman-Kittel-Kasuya-Yoshida (RKKY) interaction on nonmag- 
metic disorder and gate voltage in graphene. First the semiclassical method is employed to rederive 
the expression for RKKY interaction in clean graphene. Due to the pseudogap at Dirac point, the 
RKKY coupling in undoped graphene is found to be proportional to Next, we investigate 

how the RKKY interaction depends on nonmagnetic disorder strength and gate voltage by studying 
numerically the Anderson tight-binding model on a honeycomb lattice. We observe that the RKKY 
interaction along the armchair direction is more robust to nonmagnetic disorder than in other direc- 
tions. This effect can be explained semiclassically: The presence of multiple shortest paths between 
two lattice sites in the armchair directions is found to be responsible for the reduced disorder sen- 
sitivity. We also present the distribution of the RKKY interaction for the zigzag and armchair 
directions. We identify three different shapes of the distributions which are repeated periodically 
along the zigzag direction, while only one kind, and more narrow distribution, is observed along the 
armchair direction. Moreover, we find that the distribution of amplitudes of the RKKY interaction 
crosses over from a non- Gaussian shape with very long tails to a completely log-normal distribution 
when increasing the nonmagnetic disorder strength. The width of the log-normal distribution is 
found to linearly increase with the strength of disorder, in agreement with analytical predictions. 
At finite gate voltage near the Dirac point, Friedel oscillation appears in addition to the oscillation 
from the interference between two Dirac points. This results in a beating pattern. We study how 
these beating patterns are effected by the nonmagnetic disorder in doped graphene. 



I. INTRODUCTION 



A large effort has been devoted to the study the elec- 
tronic transport properties of graphene due to the un- 
usual nature of the quasiparticles in this material, which 
are massless chiral Dirac fermions. A recent experiment 
indicating that vacancies in graphene may induce local 
magnetic moments^ renewed the interest in investigating 
magnetic properties as well. It is belived that the carrier- 
mediated effective exchange coupling between local mo- 
ments, or Ruderman-Kittel-Kasuya-Yoshida (RKKY) in- 
teraction may play a crucial role in establishing how mag- 
netism develops in graphene. Probing these properties lo- 
cally is not completely out of reach: Using spin-polarized 
scanning tunnehng spectroscopy, a direct measurement of 
the RKKY interaction in a conventional metal has been 
done by measuring the magnetization curves of individual 
magnetic atoms adsorbed on a platinum surface.^ Similar 
experiments may soon be possible on graphene. 

Analytical and numerical studies of the RKKY inter- 
action in clean graphene at the neutrality point have 
been reported In this context, a dominant feature of 
graphene is the bipartite nature of its honeycomb lattice. 
Due to particle-hole symmetry of bipartite lattices, the 
RKKY interaction always induces ferromagnetic correla- 
tion between the magnetic impurities on the same sublat- 
tice, but antiferromagnetic correlation between the ones 
on different sublattices.^ At the neutrality point, the de- 
pendence of the RKKY interaction on the distance R be- 



tween two local magnetic moments in graphene is found 
to be in contrast to the standard behavior in con- 

ventional two-dimensional metallic systems where 1/R'^ 
is expected. In doped graphene, but not too far from 
the neutrality point, the behavior changes to 1/R^ and 
two different length scales control the RKKY interaction: 
The wavelength corresponding to the inverse of the dis- 
tance between the two Dirac points in momentum space, 
\K — and the Fermi wavelength kp^.^ 

Since the RKKY interaction is mediated by itinerant 
electrons in the host metal, nonmagnetic defects influ- 
ence directly these interactions. On-site potential fluctu- 
ations scatter the phase of the electron's wave function 
as well as cause random changes in its amplitude, alter- 
ing any interference phenomenon observed in the clean 
system. The effects of weak nonmagnetic disorder on the 
exchange interactions in conventional metals have been 
carefully studied Numerical work in the strongly lo- 
calized regime has been done in one-dimensional disor- 
dered chains These studies found that in the diffusive 
regime the RKKY interaction, when averaged over dis- 
order configurations, decays exponentially over distances 
exceeding the mean free path scale /g- However, its geo- 
metrical average value has the same amplitude as in the 
clean limit. It is only in the localized regime, on length 
scales ^ ^, i.e., exceeding the localization length ^, 
that its geometrical averaged values are exponentially 
suppressed so that the RKKY interaction is cutoff over 
distances exceeding . In our previous work we have stud- 
ied the effect of nonmagnetic disorder on the RKKY in- 
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FIG. 1. (Color online) Schematic diagrams of (a) the prop- 
agating paths of an electron in clean system and (b) the 
graphene lattice. The two sublattices are denoted as A and B 
and the two representative directions (zigzag and armchair) 
are indicated as dashed gray lines. Or is the angle between 
the displacement vector of the magnetic impurity R and the 
unit vector x. 



teraction in undoped graphene using the kernel polyno- 
mial method (KPM).^^ The unexpected and interesting 
result found was that, at weak disorder, the RKKY in- 
teraction along the armchair direction is not influenced 
by nonmagnetic disorder as much as in the other direc- 
tions. 

Motivated by this unexplained behavior, in this paper 
we start by employing the semiclassical method to evalu- 
ate the RKKY interaction in graphene. This calculation 
helped us understand the i?-dependence and the origin 
of the direction dependence of the sensitivity to disorder 
seen in Ref. 12. We then used the numerical KPM to 
calculate the RKKY interaction in graphene in order to 
study the effect of disorder and of doping in more detail. 

In order to study the disordered conduction electrons 
in graphene numerically, we employ the single-band An- 
derson tight-binding model on a honeycomb lattice. 



H = -t^^ c\cj + ^^{wi - fi) clci, 



(1) 



id) 



where t is the hopping integral, q is an annihilation op- 
erator of an electron at sitei, c] is the corresponding cre- 
ation operator, Wi is the uncorrelated on-site disorder po- 
tential distributed uniformly between [—1^/2,1^/2], and 
(z, j) indicates nearest-neighbor sites. Periodic boundary 
conditions are used and the lattice constant a and h are 
set to be unity in all calculations. 

Below, we start with the derivation of the RKKY in- 
teraction expression in clean graphene using the semi- 
classical method. 



II. SEMICLASSICAL APPROACH TO THE 
RKKY INTERACTION 

Bergmann interpreted the RKKY oscillation as an 
interference of the conduction electron's wave function 
scattered by the magnetic impurity and calculated the 
interference using the Huy gens' principle in a three- 
dimensional metaL^ According to the Huygens' principle 



in two dimensions^ the amplitude of a wave which prop- 
agates from a source at position R decays with distance 
and gains a phase factor at a position R given by 



^ik-{R-R') 



(2) 



where "^{R') is the amplitude of the wave at the source 
and the extra factor comes from the asymptotic form of 
the Bessel function in two dimensions (e*^^/y^). If an 
electron propagates from R to the origin in graphene 
(Fig.[T^), the amplitude gets a phase factor 



A = 



'^i{K+q)-R _^ ^i{K'+q)-R 



(3) 



where the wave vector is expanded around the two neigh- 
boring Dirac points K and K' with a relative wave vec- 
tor q. During its propagation, the electron gets another 
phase factor, e*^*?^/^, where Sq = vpq is the kinetic en- 
ergy of the electron near the Dirac point, 'Ui? is the Fermi 
velocity and t/2 is the propagation time from R to the 
origin. After the electron is scattered by a magnetic im- 
purity at the origin, it travels back to the position R. 
The amplitude then gains an additional modulation of 
the amplitude which is given by 



B = 



-iiK+q)-R 



i{K'+q)-R 



(4) 



where 5q depends on the properties of the magnetic im- 
purity at the origin. When the electron goes back to R^ 
its direction is opposite compared to the first propagation 
and this is the reason why the signs are opposite in the 
phase factors in A and B. In the time domain, however, 
it propagates in the same direction so that the phase fac- 
tor related with the time {eqt/2) has the same sign. This 
is consistent with the diagrammatic expression for the 
RKKY interaction which has two retarded Green's func- 
tions. For a closed path, a loop of the opposite direction 
makes a contribution with the same weight, so that the 
modulation of the charge density of an electron which 
has the energy Sq is given by 



Ap^(i^) = 2\^{R)\^AB 

_ 4(5o l + cos[(K-KO■i^] 
V R 
such that the total charge modulation is given by 

p{R) = j def{e)N{e)Ape{R) 

_ Sovj, 1 + cos[{K - K') ■ R] 



(5) 



V 



(6) 



Here, |?/^(i^)p = 1/V for free electron, V is the volume of 
the system, f{s) is the Fermi-Dirac distribution function, 
N{e) = l^l" is the density of state at the Fermi level, a is 
the pseudogap exponent, d is the spatial dimension, and 
t = 2r/vF is the total time it takes to return to R. 
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The total charge density in (Eq. (j6])) can be easily re- 
lated to the RKKY interaction.— Since in graphene d = 2 
and a = 1, the resulting distance dependence is l/R^, 
which is consistent with previous works. We conclude 
therefore that the existence of the pseudogap at the Dirac 
point causes the faster spatial decay of the RKKY in- 
teractions in graphene compared to conventional two- 
dimensional system = 2, a = 0). Even though Eq. (|6j) 
does not give detailed information of about the RKKY 
interactions, such as the sign dependence on the sublat- 
tice and an extra phase factor depending on the direction 
between point R and the origin, it agrees with previ- 
ous results obtained by the Green's function method^i^ 
namely. 



Jaa — — J 



jO _ j: 



2 1 + cos[(K - K') ■ R] 

2 3 + 3 cos[(K - • i^ + TT - 2l9i?] 
^3 ' 



(7) 
(8) 



where Or is an angle between the position vector R and 
the zigzag- AA direction {6r = 0) as described in Fig.[T]3. 



III. NUMERICAL METHOD 

We start with a general expression for the RKKY ex- 
change coupling constant in terms of the off-diagonal 
matrix elements of the local density of states pij{E) = 
{i\5{E-H)\j)^ 



-^RKKY = - J 



1) 



252 



JE<ij. J E' 



dE 



> FiE,E') 
E-E' ' 



(9) 

where F{E,E') = Re[pij{E)pji{E^)], J is the local cou- 
pling constant between the localized magnetic impurity 
and the itinerant electrons in the host metal, S is the 
magnitude of the impurity spin, i^j are the site indices 
of magnetic impurities located at position Ri (Rj)^ and 
/i is the Fermi energy. Zero temperature (T = 0) is as- 
sumed. Using the KPM, one may calculate the matrix 
elements of the local density of state efficiently^ 



1 



Pij 



ttVI - E'^ 



M 



go rn^ 



2^gimi^Ti{E) 



1=1 



(10) 



where Ti{E) is the Ith Chebyshev polynomial of the first 
kind, 171^ = {i\Ti{H)\j)^ H is an electronic Hamiltonian, 
and gi are the Jackson kernels coefficients. The sum is 
taken up to a cutoff number M. One may obtain the 
expansion coefficients my using the recurrence relation 
of Chebyshev polynomials, namely, T^+i = 2HTi{H) — 
Ti_i{H) with To{H) = 1 and Ti{H) = H. Implicit in 
Eq. (p^Oj) is the normalization of the energy spectrum to 
a band of unity width. 



IV. DISORDERED GRAPHENE AT HALF 
FILLING (/i = 0, / 0) 

We have investigated the effect of on-site nonmagnetic 
disorder at the charge neutrality point {p = 0^ W ^ 0). 
For each value of 1600 different disorder configura- 
tions were generated and the corresponding matrix el- 
ements of the density of state pij (Eq.[TQ|) evaluated 
through the KPM with a Chebyshev polynomial cut- 
off number M = 5 x 10^ on a lattice with 5 x 10^ 
sites. We have previously reported that in the diffusive 
regime (/e < R < which was characterized by deter- 
mining the mean free path /g and the localization length 
^, the armchair direction is the only direction in which 
the averaged value of the RKKY interaction over disorder 
configurations does not alter its sign.^^ To illustrate this 
effect, we calculated the RKKY interactions along the 
directions {Or = 7r/12.5, tt/IO) which pass through only 
the same sublattice and the results are shown in Fig. [2] 
together with the ones for the zigzag- AA and armchair- 
AA directions. Notice that the amplitudes of the RKKY 
interaction are multiplied by the cube of the distance in 
order to make the effect of nonmagnetic disorder more 
transparent. As expected from previous studies^ — the 
amplitude of the averaged interactions decays exponen- 
tially over length scales exceeding the mean free path /g. 



(-^RKKY)avg = -^RKKY ^ 



-R/le 



(11) 



where (O)avg indicates averaging over disorder configura- 
tions and Jrxky amplitude of the RKKY interac- 
tion in a clean system. On-site disorder breaks the sub- 
lattice symmetry. Consequently, the sign of the RKKY 
interaction between the moments which are localized at 
the same sublattice fluctuates, allowing both ferromag- 
netic and antiferromagnetic correlation (Fig.[2^, b, c). 
The importance of the sublattice symmetry was high- 
lighted in our previous work by considering the hopping 
disorder When randomness is added to the hopping 
integral t = to -\- At without on-site disorder, it does 
not break the sublattice symmetry and the RKKY inter- 
action never changes its sign for magnetic moments sit- 
ting in either the same or different sublattices, even for a 
fixed disorder configuration. Interestingly, the averaged 
RKKY interaction amplitude along the armchair direc- 
tion with on-site disorder does not change sign (Fig.[2]i), 
while for a particular disorder configuration it randomly 
changes both sign and amplitude. 

In a dirty metal the charge density is modified by dis- 
order due to the random phase factors. When these ran- 
dom phase factors are simply averaged over disorder con- 
figurations, they gain an exponentially decaying factor. 
Due to the geometrical anisotropy in graphene, the in- 
fluence of these random phase factors depends on the 
path direction. According to Feynman's path integral 
representation, the coupling is dominated by the short- 
est path between two magnetic impurities. There always 
exists an even number of shortest paths along the arm- 
chair direction (Fig.[3]3), but only one along the zigzag 
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FIG. 3. (Color online) Schematic diagrams of the short- 
est path from the origin to the first lattice point along (a) 
the zigzag-AA direction and (b) the armchair- AA direction. 
^1 (2) denotes the phase shift due to nonmagnetic disorder and 
A (B) denote the modulation of the charge density amplitude. 



FIG. 2. (Color online) Plots of the RKKY interaction 
strength multiplied by the cube of the distance, R^, along 
the (a) Or = (zigzag-AA) (b) Or = 7r/12.5 (c) Or = tt/IO 
and (d) Or = tt/S (armchair- AA) direction in the diffusive 
regime, as averaged over 1600 different disorder configura- 
tions. A lattice with 5 x 10^ sites and a polynomial degree 
cutoff of M = 5 X 10"^ are used in these numerical calculations. 



direction (Fig.[3ti). Therefore, for the zigzag direction, 
the electron wave is scattered by the same disorder twice 
with the same phase (e.g. Si in Fig. [3^) when it returns 
to the origin. However, along the armchair direction, 
there are closed paths which include different impurities 
and thereby different scattering phases (e.g. and 62 
in Fig.inb). If we average the modulation of the charge 
density over the phases Si and 62 we find 



-2{5i) 



and 



-2(5l) 



(12) 



(13) 



where (• • •) is the average over disorder configurations 
and Apo is the modulation of the charge density in 
the clean system [Eq. (|6|)]. The modulation along the 
armchair direction [Eq. (p!3|) ] is dominated by the term 
Apo e~^^^^ coming from the closed loop connected by two 
unrelated paths, so that the average value is exponen- 
tially closer to the clean value, while the zigzag direction 
is dominated by Apo e~'^^^^^ [Eq. (p!2]) ]. In this sense, one 
may expect the RKKY interaction in the armchair direc- 
tion to be less sensitive to the nonmagnetic weak disorder 
than in other directions. This may explain our previous 
result^ and the results shown in Fig.O where the av- 
eraged value in the armchair direction remains susbstan- 
tially larger than in the other directions and does not 
change its sign. 

In order to investigate in more detail how much non- 
magnetic disorder affects the RKKY interaction and how 
the latter depends on direction in the graphene lattice 
and on particular lattice points, we have calculated the 



distribution of the interaction amplitude for a disorder 
strength W = t and at lattice points which are on the 
zigzag or armchair direction. The results are shown in 
Fig.|4]a-e. In these calculations, 3 x 10^ different disorder 
configurations are used with for a lattice of 2 x 10^ sites, 
while the number of Chebyshev polynomials M = 10^ is 
fixed. When the distance R between the magnetic im- 
purities is smaller than the mean free path, the RKKY 
interaction multiplied by the cube of the distance has a 
shape which does not depend on the distance. Interest- 
ingly, the distribution has three different shapes along 
the zigzag-AA direction (Fig.|4^-c), which are repeated 
periodically every third site in that direction. However, 
there is only one type of distribution along the armchair 
direction, as can be seen in (Fig.SJi). The oscillating 
factor in the semiclassical expressions for the charge den- 
sity, cos[(i^ — K^) • R] in Eq. (|6]), takes only the values 
,(1,-1/2,-1/2) along sites in the zigzag-AA direction. 
The RKKY interaction on sites which give the same nu- 
merical factor, either 1 or —1/2, are found to have the 
same distribution. In order to directly compare them 
with each other, we plot together the distributions of the 
RKKY interactions obtained for these sites in Fig.H^. 
The RKKY interaction at sites where the oscillating fac- 
tor takes the value (-1/2) along the zigzag-AA direction 
has the broadest (green) distribution, while the RKKY 
interaction at site where the oscillating factor takes the 
value (1) has the narrowest (red) distribution. 

The unaveraged RKKY interaction, i.e., that obtained 
for a particular disorder representation, is found to re- 
main long ranged as in the clean system. However, we 
find that the oscillations acquire a random phase that 
also changes with the distance R between magnetic impu- 
rities. Therefore, the average value does not characterize 
the RKKY interaction at distances exceeding the mean 
free path. To find the typical value, we evaluated the 
geometrical average of the interaction amplitude, which 
is defined by 



xgeo 
^RKKY 



exp 



^(^^(JrkKy)^ )avg 



(14) 



The typical value is found to have the same power- 
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FIG. 4. (Color online) Plots of the distribution of the RKKY 
interaction amplitude ( ^rkky) multiplied by the cube of 
distance, i?^, for the (a) first (b) second (c) third of the triplet 
of points along zigzag- AA direction (see Fig. 4 f for the def- 
inition of these points), and (d) for the armchair direction. 
The disorder strength is fixed to VK = t. For comparison, 
these distributions are plotted together in (e). A lattice with 
2 X 10^ sites and a polynomial degree cutoff of M = 10^ are 
used. The lattice constant a is set to unity. 



FIG. 5. (Color online) (a) The averaged density of states 
in a double logarithmic plot (b) and the geometrical aver- 
aged RKKY interactions along the zigzag- A A direction. For 
the density of state calculation, a lattice with 3 x 10^ sites 
and a polynomial degree cutoff of M = 7 x 10"^ are used. 
For the RKKY interaction calculation, the same lattice size 
and polynomial cutoff of Fig. [2] are used. The black lines 
in (a) represent fittings to the relation p{E) — j\E\ with 
7 = 0.95, 1.07, 1.3 for Vl^ = 0, 0.5t, t, respectively. 



law decaying behavior with distance as the amplitude of 
the RKKY interaction in the clean limit A"— i^^ii^ Before 
a direct evaluation of the geometrical average, we have 
calculated the density of state for two disorder strengths 
W = 0.5t, t to observe how weak disorder affects the 
pseudogap at the neutrality point {E = 0), since, as 
we have seen above, the power of the pseudogap is di- 
rectly related to the unconventional power-law decay of 
the RKKY interaction. For the density-of-states calcu- 
lation, we used the KPMi^iiS method with 3 x lO'^ sites 
and a polynomial degree cutoff of M = 7 x 10^. As one 
may see from the plot of the density of states in (Fig. [5^), 
the pseudogap is still not filled for weak disorder, has the 
same power law as in the clean limit, and only the slope 
around the neutrality point is changed, 

P{E) = 1\E\, (15) 

where p{E) is the density of states and E is the energy 
measured in units of the hopping amplitude t. The slope 
7 depends on the disorder strength and is obtained by 
fitting the data in Fig. [6l 

Using the Born and T-matrix approximations, an an- 
alytical study has reported that there is a logarithmic 



correction to the density of states around the neutral- 
ity point, yielding p{E) = |£^| In |£^| in the presence of 
disorder, This logarithm correction does not change 
the power of the distance dependence of the RKKY in- 
teraction [Eq.®]. Consequently, the geometrical average 
of the RKKY interaction is expected to have the same 
power-law exponent as the clean system (1/i?^). The di- 
rect numerical calculation shown in Fig.lSb strongly sup- 
ports this conclusion. 

Figure[6^ shows how the RKKY interaction in a 
strongly disordered sample gets suppressed by several 
orders of magnitude when the strength of the nonmag- 
netic random potential is increased. In order to inves- 
tigate the broadening of the distribution, we employed 
3 X 10^ realizations of the disorder potential and calcu- 
lated the RKKY interaction amplitude for a fixed dis- 
tance R = 50V^. The results are shown in Fig.[6]3. 
The inset indicates that the interaction strength fol- 
lows a distribution with very long assymetric tails. The 
squared amplitude (-^rkky) ^ distribution similar to 
log-normal with a width that increases with disorder 
strength. Using a field-theoretical approach valid in the 
metallic regime, Lerner found^^ that the increase in the 
strength of the nonmagnetic disorder leads to a crossover 
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FIG. 6. (Color online) Plots of (a) the absolute value of the 
RKKY interaction in a disordered sample as a function of 
distance and (b) of the distribution of the logarithm of the 
RKKY interaction amplitude for R — 50A/3tt and different 
disorder strengths. The inset of (b) shows the distribution 
of the amplitude itself for W = 41. A lattice with 1.8 x 10^ 
sites and a polynomial degree cutoff of M = 3 x 10^ are used. 
The black dashed lines represent the resulting distribution 
functions with a = 2.4, 3.2, 4.35 and a = 32.5, 50, 72.5 for 
W = 4t, 6t, 8t, respectively (see text). 



in the shape of the distribution function from a broad 
non- Gaussian with the very long tails in the weak disor- 
der regime to a completely log-normal distribution in the 
region of strong disorder regime. This may explain the 
crossover of the distribution which we observe from weak 
disorder W = t to strong disorder {W = 4t, 6t, St). In 
order to directly compare we have fitted the results with 
a log-normal functional form, which is given by 



P{x) 



N 



exp 



{x — 



2a2 



(16) 



where TV = 3 x 10^ is the number of realizations and 
X = ^^[J^kky]/"^' data are shown together with the 
fitting curves in Fig.[6]3 (black dashed lines). 

The width cr of the distribution in Eq. (p!6|) has been 
analyzed as functions of the disorder strength W and is 
shown in Fig. [71 The red line in Fig. [3 is a linear fitting 
curve (cr = W/2 + 0.34) and agrees with the analytical 
prediction by Lerner, who has used the renormalization 
group method to obtain^^ 



cr ~ 1 / ^ W, 




FIG. 7. (Color online) Plot of the width of the distribution 
of RKKY interaction amplitudes as a function of the disorder 
strength (in unit of t). The red line represents a linear 
fitting curve which yields a — W/2 + 0.34. 



where le ~ is the mean free path of electrons, 

which we studied in our previous workJ^ 



V. DOPED GRAPHENE (/i / 0) 

A. Clean system (W = 0) 

By controlling /i in the Hamiltonian of Eq. ([T]) , we in- 
vestigate how the RKKY interaction evolves with the 
Fermi level. Results are shown in Fig.[8l where the in- 
teraction amplitude is multiplied by i?^ in order to em- 
phasize its oscillatory behavior. In these calculations, we 
used a lattice with 7.2 x 10^ sites and a polynomial degree 
cutoff M = 3 X 10^. Near the neutrality point a beating 
pattern appears as shown in Fig. [8^, b. It consists of a 
superposition of waves with the wave vectors K — K' and 
qf^ where gr^? is the Fermi wave vector originating from 
the Friedel oscillations at finite Fermi energy. Recently, 
the following analytical expressions for the beating pat- 
tern were derived using lattice Green's functions:^ 



AA 



and 



-^AB = -^A 



AB 



SqpR 



37r 



-2,1 
2,4 



1 3 

2 ' 2 

1,2,0,- 



(18) 



, (19) 



(17) 



where ^aa(b) RKKY coupling function at the neu- 

trality point [see Eq. (j8j)] and G is the Meijer-G function. 
Note that term within brackets describes the isotropic 
dependence of the oscillations on the Fermi momen- 
tum qp. The external prefactor, -^aa(B)' contrary, 
is strongly anisotropic, depending on the vector given by 
the momentum difference between the two neighbored 
Dirac points K — K' . To make a comparison with our 
calculations, the function represented by Eq. (p!9|) is also 
presented in Fig. [8] by a red dashed line. Excellent agree- 
ment is found. One can estimate the wavelength of the 
long oscillation, which appears at finite Fermi level, using 
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FIG. 8. (Color online) (a) The RKKY interaction amplitude 
multiplied by in doped graphene at the chemical potential 
fi — O.lt along the zigzag- AB direction. Density plots of the 
RKKY interaction amplitude for (b) /i = 0.2t and (c) /i = t 
for all directions. In (a), the result of culations used the kernel 
polynomial method and the lattice Green's function method 
are represented as solid blue and dashed red line, respectively. 
A lattice with 7.2 x 10^ sites and a polynomial degree cutoff 
of M = 5 X 10^ were used. The lattice constant a is set to 
unity. 



the dispersion relation near the neutrality point, which 
is given by 



E{q) = VF\q\ 



(20) 



where q is the momentum relative to the Dirac point 
and vf = 3ta/2 is the Fermi velocity ^21 The Fermi wave- 
length is found to be about ^ 50a for /i = O.lt. This 
coincides with the period seen in (Fig. [8^). As expected 
from Eq. ([19]), the oscillations with large period seen in 
(Fig.[8]3) are isotropic. 

We have also calculated the RKKY interaction am- 
plitude for highly doped graphene (/i = t). The results 
are multiplied by the square of the distance, i?^, and 
are shown in the density plots of (Fig.[8t). The behav- 
ior cannot be described by Eqs. (fTSj) and ([19]) which are 
valid only close to the neutrality point. When the Fermi 
level is exactly at the van Hove singularity at (/i = t), 
the ordering pattern of the RKKY interactions along the 
zigzag direction is reversed. In other words, the correla- 
tion between impurities on zigzag- AA or BB is always an- 
tiferromagnetic and ferromagnetic for zigzag- AB or BA 
pairs. At the same time, the interactions are strongly 
suppressed for the other directions. This is in accordance 
with a result of a previous study 
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FIG. 9. (Color online) The average RKKY interaction 
strength multiplied by the square of the distance, i?^, along 
the (a) zigzag-AA and (b) armchair-AA directions at the 
Fermi energy /i = 0.2t is plotted in the diffusive regime. 1600 
different disorder configurations are used in the averaging. 
The same lattice size and polynomial cutoff as in Fig. [2] were 
used. 



B. Disordered system (W 7^ 0) 

We have also calculated the RKKY interaction ampli- 
tude in disordered doped graphene and the results are 
shown in Fig.|9l A lattice with 5 x 10^ sites and a poly- 
nomial degree cutoff of M = 5 x 10^ were used in these 
numerical calculations. In our previous work^ we re- 
ported that the amplitude of the averaged RKKY in- 
teraction along the zigzag-AA direction in the ballistic 
regime (i? < /g) increases with weak disorder (Fig. |4]F). 
Even though the unusual oscillation coming from the in- 
terference of two Dirac points still exists (Fig. (9^), the 
amplitude of the envelope of the averaged RKKY inter- 
action decreases with disorder strength W and the pe- 
riod (27r//cF) of the oscillation coming from the finite 
Fermi energy is modified by the disorder as in a con- 
ventional metal. 



VI. DISCUSSION AND CONCLUSIONS 

In conclusion, we have employed a semiclassical de- 
scription of the RKKY interaction in terms of the modu- 
lation of the charge density by the presence of magnetic 
impurities in order to rederive an analytical expression 
applicable to pure graphene. This semiclassical approach 
provides not only a simpler derivation, but also an intu- 
itive interpretation of the origin of the oscillating behav- 
ior of the RKKY interaction at zero doping as an interfer- 
ence of the two degenerate Dirac points. Moreover, the 
origin of the unusual power-law decay of the RKKY inter- 
action in pure graphene at the neutrality point is clearly 
related to the pseudogap in the density of states. Using 
a Feynman's path integral scheme, we could also trace 
the origin of the anisotropic sensitivity of the RKKY in- 
teraction to disorder to the presence of multiple shortest 
paths between two magnetic impurities in the armchair 
direction, showing that this direction is more robust to 
disorder than the zigzag one. 

As an extension of our previous study^ we have cal- 
culated and studied the RKKY interaction in doped and 
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disordered graphene in detail using the kernel polynomial 
method. Fnite gate voltage breaks particle-hole symme- 
try and the resulting finite Fermi surface yields Friedel 
oscillations, so that the sign of the RKKY interaction be- 
tween the impurities localized on the same sublattice now 
oscillates with distance. When the Fermi level is exactly 
at the van Hove singularity (/i = t), the ordering pattern 
of the RKKY interactions along the zigzag direction is 
reversed. At the same time, the interactions are strongly 
suppressed for all other directions. 

In order to study the anisotropic infiuence of non- 
magnetic disorder, we evaluated the RKKY interactions 
along two different directions between the zigzag and 
armchair directions. As reported in our previous study 
and expected from the semiclassical approach, in the dif- 
fusive regime the armchair direction is not affected by 
the nonmagnetic disorder as much as the other direc- 
tions. We have also found that, in the ballistic regime 
{R < le)^ the distribution of the RKKY interactions 
along the zigzag direction is not universal but depends on 
the lattice sites at which the pair of magnetic impurities 
sit. We identified three different representative shapes, 
which repeat themselves periodically. By an accurate 
evaluation of the density of states around the neutrality 
point in weakly disordered regime {W < t), we confirmed 
that the linear dispersion relation is still valid and the 
pseudogap is not filled. This is in full agreement with 
the fact that the geometrical average of the RKKY in- 
teraction in the diffusive regime decays as in the clean 
system, namely, as and not as which is the 

usual behavior for two-dimensional metals. In the lo- 
calized regime (i? > ^), the geometrical average is also 
exponentially suppressed at distances exceeding the lo- 
calization length ^ and the distribution of the strength 
of the RKKY interaction shows a crossover from the non- 
Gaussian shape with very long tails to the completely log- 
normal form when increasing the disorder strength. We 
have analyzed the width of the log-normal distribution 
and confirmed that it increases linearly with the ampli- 
tude of the disorder potential W. 

In this work we showed that the KPM method is ef- 
ficient and accurate for studying interactions in disor- 
dered systems. In order to minimize the computation 
time while keeping the highest accuracy, the convergence 
of the calculations with respect to the Chebyshev poly- 
nomial cutoff degree M and the system size L have been 
investigated in detail (see Appendix A). The proper cut- 
off M to reach convergence is found to increase linearly 
with the distance R between two magnetic impurities, as 
seen in Fig. [10] b. For a given distance R^ a system size L 
about five times larger than R is found to give convergent 
results. In comparison with the exact numerical diago- 
nalization, the KPM is found to be much faster. It can 
be implemented for very large system sizes, even in dis- 
ordered ones, where thousands of realizations are needed 
in order to yield meaningfull statistics. 
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FIG. 10. (Color online) Plots of (a) the RKKY interaction as 
function of the polynomial degree M in a lattice with 5 x 10^ 
sites, (b) The smallest cutoff number M that yields a conver- 
gent result as function of the distance between the magnetic 
impurities R. 



Appendix A: Convergence of the Kernel Polynomial 
Method Calculations 



When using the KPM, the calculation of the Cheby- 
shev polynomials using recurrence relations consumes 
most of the computation time. Therefore, we investi- 
gated first the relation between the cutoff number M 
and the convergence of the results in order to be able 
to minimize M and optimize the calculations. For clean 
graphene, a lattice with 5 x 10^ sites was used in these cal- 
culations. When a cutoff number M is not sufficient, the 
amplitude of the RKKY interaction deviates from the ex- 
pected power-law behavior, as indicated by the blue and 
red lines in Fig.fTOk). When we determined the smallest 
cutoff number M such that the variance of the ampli- 
tude of the RKKY interaction is less than 5%, we found 
that it increases linearly with the distance between two 
magnetic impurities R as seen in Fig.fTOb). This linear re- 
lation between the distance and the cutoff number allows 
the rapid calculation of the RKKY interaction amplitude 
between the magnetic moments even at large distances 
R. 

In order to minimize the KPM calculation time further, 
we have also studied the smallest system size L which 
yields convergent results, as shown in Fig.[TTJ Rp denotes 
the longest distance used in the calculation, Rp = SOV^a 
and the cutoff number is M = 5 x 10^. One can observe 
a good convergent behavior (see green and black line in 
Fig.fTT]) when the system size L is larger than bRp. The 
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FIG. 11. (Color online) Plots of the RKKY interaction in 
terms of a system sizeL. 50V3a is used as the longest dis- 



tance Rf- a cutoff number M 
calculations. 



5 X 10 is used in these 



exact diagonalization method also yields a proper result 
when L = 5Rf'^ However, the KPM does not require ma- 



trix diagonalization and therefore is a much faster com- 
putational tool. 
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